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1. Introduction 

The nonlinear Sclirodinger (NLS) equation is a ubiquitous nonlinear wave equation with 
a range of applications including the propagation of hght within a waveguide lU [2], 
the behaviour of deep water waves [3j, and the mean-field theory of Bose-Einstein 
condensates [1]. However, in many practical situations the NLS represents only the 
zeroth order approximation to the system, and for this reason, the response of an NLS 
system to small perturbations is important |5]. A novel mathematical formalism (based 
on the four linearly independent Goldstone modes of the linearised problem) with which 
one can deal with the spatial consequences of such perturbations is the aim of this 
article. Currently relevant examples of such perturbative mechanisms include the loss 
and/or dephasing of coherent light traveling through an optical fibre, and the presence 
of quantum and / or thermal noise in Bose-Einstein condensates [6] . 

The problem under consideration has received a significant amount of attention in 
previous literature [TJ El [9l [TOl ttH US Il3], and in fact it would seem that a general 
approach toward such problems has been established within the community since the 
late 1970's. Briefly, the approach focuses on finding eigenfunctions of a differential 
operator which is found from linearising the NLS around an analytic soliton solution. 
The majority of the earlier work was concerned more with the bright solitons found in the 
self- focusing NLS [H]. Progress on the dark soliton of the self-defocusing NLS caught up 
with its bright counterpart in the mid-to-late 1990 's, with the introduction of a complete 
set of so-called "squared Jost solutions" [15]. The crux for the dark soliton solutions 
involved dealing with the nonvanishing boundary conditions, this issue was avoided in 
the earlier work of Ref. [12], where they force a vanishing boundary condition onto 
the perturbation for theoretical convenience. Reference [13] develops a method based 
on separating out the internal soliton dynamics from that of the boundary conditions, 
however such a separation is approximate at best [15], HE] . These squared Jost solutions 
of Ref. [15] elegantly provided the desired eigenfunctions for all real eigenvalues, except 
the case where the eigenvalue is zero (in this case the eigenfunctions are commonly 
referred to as the Goldstone modes). In this limit as the eigenvalue tends toward zero, 
the squared Jost solutions collapse down to just two linearly independent solutions 
(the linearised differential operator is ultimately a fourth order differential equation 
and should therefore yield four linearly independent solutions). This fact was noted 
in Refs. [151 IE]; cind two additional generalised eigenvectors were introduced to cope 
with the absence of the remaining two solutions. With the inclusion of these generalised 
eigenvectors it was shown that one had a complete set of functions. The main results 
of Ref. [IS] has led onto several other publications [HI [191 EQ] EH [22] of a similar vein. 

The issue has recently seen an influx of interest coming from the community of 
scientists involved with ultra-cold quantum gases. The original observation of dark- 
soliton-excitations of Bose-Einstein condensates within elongated trapping geometries 
came in 1999 [23j and continues to accrue an impressive number of citations. 
Sophisticated numerical techniques have been employed in Refs. [211 |25l [26] which 
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investigate the lifetimes of dark solitons in the presence of quantum and thermal noise 
(special attention was paid to the high temperature regime in Ref. ^7\). Analytic 
approaches toward the same problem were put forth earlier in Refs. [2H1 EH], and 
included the effects of the anomolous modes associated with phase diffusion of the 
Bose-Einstein condensate (originally considered in Ref. [30]) as well as diffusion in 
the position of the dark soliton (these anomolous modes are given in equations UM - 
[15] of this article). Conflicting interpretations of the ensemble density evolution sparked 
debate as to whether the soliton exhibits decay or diffusion in the presence of noise |31j . 
Another mechanism put forth as being responsible for the decay of dark solitons is 
the effective three-body contact interactions considered in [32l [33]. They claim the 
soliton is protected against decay by the integrability in the system under two-body 
collisions. This integrability must be broken to observe soliton decay, a hypothesis which 
is supported by the claims of [28l [29] . The inclusion of three-body interactions destroys 
the integrability in the system. Further experiments in the field have successfully verified 
much of the fundamental interest in solitons such as their particle-like properties and 
mutual transparency under collision |3ll [35] . 

The overall goal of the present paper differs slightly from much of the previous 
literature, specifically we emphasise that time evolution of the soliton parameters is not 
addressed in this article (see Ref. [151 120] for a treatment of this problem). Rather 
we concern ourselves solely with the first order correction to the spatial profile of the 
soliton. This correction is found by solving a nonhomogeneous fourth order differential 
equation (see equation [7] in section [2] of this paper). It is true that this correction can 
in principle be dealt with using the complete set of Ref. [15], however this can be very 
difficult in general. Indeed it is expressly stated in Ref. [15] (see the final paragraph of 
the introduction) that the first order correction is difficult to obtain via their method. 
We present here a much simpler method based on analytic solutions for all four linearly 
independent Goldstone modes. It is the introduction of these analytic expressions for 
the two, previously unpublished, Goldstone modes which allows us to proceed in this 
way. The four solutions are related to the four fundamental symmetries of the NLS. 
These symmetries are phase symmetry, translational symmetry, Galilean symmetry, 
and dilaton symmetry. Aside from the methods aesthetic appeal, Ref. |36j describes a 
physical system (which necessitated the authors interest in this field), where, due to the 
numerical nature of the perturbing function, (denoted g{x) in equation [7] of the current 
paper, but denoted f{z) in Ref. [36]) the method of Ref. p3] was rendered useless. 

The paper is organised as follows: In section [2] we set up the problem by linearising 
the NLS equation around a dark soliton. In section 13.11 we look at the squared Jost 
solutions of Ref. [15] and discuss their importance as eigenfunctions of the linearised 
problem. In section [3^ we look at how these squared Jost solutions behave in the limit 
as the eigenvalue tends to zero. After establishing the fact that (in this zero eigenvalue 
limit) the squared Jost solutions give only two of the four possible eigenvectors, we 
give exact analytic solutions for all four eigenvectors. In section [H we use these 
eigenvectors to construct a Greens matrix for the differential operator of the linearised 
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problem. In section 15.11 we illustrate the use of this Greens matrix in solving a practical 
example (specifically the correction to the dark-soliton wavefunction of a Bose-Einstein 
condensate, in the presence of fluctuations). 

2. Basic formalism 

The usual nonlinear Schrodinger equation (with a defocusing nonlinearity) , in its 
dimensionless form, is 

_,5^^_l52^+|^|2^ = 0, (1) 

which, after a Galilean boost of the coordinates, {x = z — vt) becomes 

-idti)-]^dl^lj + ivd^^lj + \^lj\'^ij = (2) 

An interesting solution to equation [2] under non-vanishing boundary conditions is 
Tsuzuki's single soliton solution [371[38]. In this case, the function can be separated into 
the product ip{x,t) = e~^^ipQ{x), (the Galiliean shift is important for this separation). 
The solution is then 

'ipo{x) = cos{9) tanh (xc) + i sin(6'), (3) 

where v = sin(^) is the velocity of the soliton, and we have introduced a position 
coordinate — X COS (6) for notational convenience. The boundary condition in use is 
^ 1 as X ±oo (i.e. |V^oP is normalised to unity far away from the soliton). 
Now let us consider a perturbation to this NLS system of the form; 

- idt^ - ^dl^/j + ivd^i, + = eF\^,^\. (4) 

where < e ^ 1 and F\i\}^i\}\ represents some process responsible for the departure 
from the ideal NLS and denotes a complex conjugate. In a similar vein to Tsuzuki's 
solution of the unperturbed solution, we seek a separable solution in the form 

^(x,t) =e-** [V^o(a;) + e7/^i(x,To,Ti,...) + e27/'2(x,To,Ti,...)] (5) 

where the coordinates T„ = e"t, for n = 0,1,2,..., introduce a multiple-time-scale 
analysis. In the limit as e — )■ the coordinates To,Ti, . . . may be regarded as being 
independent. 

As an aside, we note that a solution to equation H] in the form of equation [5] is 
certainly not guaranteed, however the ansatz may be appropriate in certain scenarios. 
To aid any reader, who is interested in the application of this work, in determining 
whether or not the ansatz of equation [5] is appropriate in a particular case we outline a 
few basic points. 

• When e = the system is a perfect NLS system and the function ip is given by 
Tsuzuki's single soliton solution. Changes in i{j occur over a length scale Xc ~ 1 
and a time scale t ^ 1. 
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• For finite e the system will acquire an additional dynamical evolution which occurs 
over a timescale et ~ 1 |9], aswell as a new spatial profile (given by the spatial 
dependence of ipi) which is an 0(e) correction to ipo{x). 

Continuing on with the formalism, we expand the time derivative as dt = 9To+ec^Ti + 
. . . and look for a solution of ipi under the assumption that the rapid-time evolution 
(if any exists) is complete, that is dxoi'i = 0. Inserting equation [5] into equation H] and 
keeping only the terms which are linear in e we get 



--Dl + tvD, + 2\^o\ 



(6) 



where Da = Crucially for this particular approach to be relevant, the right-hand- 
side of equation [H] should not depend on the rapid-time variable Tq. The severity of this 
condition is unclear in general, however, at least in the case of one- dimensional Bose- 
Einstein condensates (where the author first encountered this kind of problem) this 
condition is certainly true. The problem then is finding a solution for the perturbation 
ipi. This is given by the following fourth order, nonhomogeneous differential equation; 



^1 










^1 






9{x, 


Ti) _ 



(7) 



where 



(8) 



The function g is the right hand side of equation [6l and can only depend on the slow- 
time variable Ti. We will refer to the linear operator Hx as the linearised operator. The 
eigenfunctions of this operator play an important part in the solution to equation [71 



3. Eigenfunctions of the linearised operator 

3.1. Non-zero eigenvalues 

In this section we briefly review some previous literature on this problem 
Specifically we look for solutions to 



ue{x) 


= E 


ue{x) 


ve{x) 




-ve{x) 



(9) 



for a fixed E ^ Four linearly independent functions and v-'e can be found by 
searching the previous literature 



Ui 



e'"^'' [kj/2 + E/kj + i cos{e) tanh (xj]^ 



(10) 
(11) 



where j = 1, 2, 3, 4 and kj is one of the four roots to the polynomial [E + k sin( 
k'^{k'^/4 +1). It is worth while to note that two of the roots {ki and ^2 say) are real, 
while two of the roots (^3 and k^ say) are complex. The complex roots mean and 
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diverge exponentially as x tends to either positive or negative infinity and for this 
reason are usually excluded on the grounds that they are unphysical. 

Equations [T014TT] can be thought of as the radiative eigenvectors of Tix- Plane wave 
excitations moving through the system essentially see the dark soliton as a reflectionless 
potential and emerge on the other side with nothing more than a phase shift. 



3.2. Zero eigenvalues 



As well as the radiative eigenvectors of the previous subsection, one also has a discrete 
set of eigenvectors associated with the symmetries of equation [H These are nonradiative 
eigenvectors and are commonly referred to as Goldstone modes. They have zero energy, 
but they have physical effects such as changing the phase of the soliton, shifting its 
spatial position, or dilating its profile. We thus turn our attention to solving the 
homogeneous problem. 



IjJ{x) 







U!{x) 








(12) 



to find these Goldstone modes. The fact that equation [9] is solved for 7^ would 
seem to indicate that solutions to equation [12] could be found simply by taking the limit 
-E — 7- 0. Unfortunately this isn't the case, as -E — > the four solutions of equations 
[TO] - [TT] collapse down into just two linearly independent solutions. 





(x) 




i {cos{6) tanh (xc) + i sin(6')) 




iipo 


U!i 


(x) 




—2 (cos(^^) tanh (xc) — i sin(6')) 






W2( 


x) 




sech^ (Xc) 








^2( 


x) _ 




sech^ (xc) 









and so we find that two of the solutions are absent from the previous literature. This 
point has not gone unnoticed, and the usual strategy for dealing with these absent 
solutions is to find generalised eigenvectors which satisfy 



'Hx 



n{x) 




uj{x) 


0(x) 




U!{x) 



(13) 



The previous literature contains expressions for two such generalised eigenvectors (see 
for example, appendix A of Ref. [17]) and it is the union of the T-L^ and Til null-spaces 
which is then used to form a complete set of functions. 

Rather than adopt this approach based on generalised eigenvectors, we write down 
expressions for all four linearly independent solutions to equation [12] 

uji{x) = — sin(^^) + i cos{9) tanh(xc) (14) 
ui2{x) = sech^(xc) (15) 
u!2,{x) = sech^(xc) [2xc — XcCosh(2xc) + (3/2) sinh(2xc)] tan(6') + 

2i [xctanh(xc) — 1] (16) 
uj^ix) = sech^(xc) {xc (10 - 4cos^(6') - 8sin(6') sin(6' - 2iXc)) + 
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cosh(a;c) [i sin(26' — 3iXc) — 5i sm{29 — ixc)] + 6 sinh(2a;c)} • (17) 

These four expressions form the key resuh of this paper {ui and uj2 have appeared in the 
previous hterature, however to the best of our knowledge and have not). These 
expressions do not follow from the finite E eigenvectors, rather they are related to the 
four fundamental symmetries of the NLS; Ui ^ phase symmetry, uj2 ^ translational 
symmetry, -H- Galilean symmetry, and UJ4 -H- dilaton symmetry. A brief summary of 
these symmetries is given below: Assuming that (f)o{x,t) is a solution of equation [1] and 
a is any real constant, then 

• phase symmetry tells us that 0o(x, t) = e'^°'(f)Q{x,t) will also be a solution, 

• translational symmetry tells us that (f)Q{x,t) = (f)o{x — a,t) will also be a solution, 

2 

• Galilean symmetry tells us that 0Q(x,t) = e^^°'^~^^^(f)o{x — at,t) will also be a 
solution, 

• dilaton symmetry tells us that (/)Q(x,t) = a(l)o{ax , a'^t) will also be a solution. 

In order to show the linear independence of equations [HHTT] we calculate the 
Wronskian 

UJi{x) U2{x) Us{x) U4{x) 

uj[{x) ^2(0;) 1^3(3;) 1^4(3;) 
io'lix) io'^ix) w^'(x) <(x) 
tu'l'ix) <(x) u'^'ix) <(x) 

and we see that, provided < 9 < 7i/2, the solutions are hnearly independent. In the 
case where 6 = tt/2 the soliton has vanished from the system and the problem becomes 
trivial. 

4. Constructing a Greens matrix 

Returning our attention to the solution of equation [3, we use the zero-eigenvalue 
solutions given in equations [T4l - [T7l to construct a Greens matrix for the linearised 
operator. The minimum requirement for this Greens matrix being that it satisfies the 
following condition; 

n^G{x,s) = l2S{x - s) (19) 

where I2 is the 2x2 identity matrix and G denotes the 2x2 Greens matrix. The general 
solution to equation [7] will then be given by 

V^i(x) 

Additional requirements given by symmetry and boundary conditions of the specific 
problem will completely determine G. 



= 512 cos^(^) sech^(a;c) sin^(^ - ix^), (18) 



G{x, s) 



(20) 
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G{x,s) = Y, 
i=i 



UJj{x) 



Kj[S) K,j[S) 



Xj{s) Xj{s) 



S < X 



X < s 



(21) 



and equation [19] gives rise to the conditions when x = s, 

hm G{x, s) = hm G{x, s) 



hm DxG{x, s) 



hm DxG{x, s) 



-2h 



These conditions manifest in the foUowing simuhaneous equations for kj and A 



'J' 



Ki{s) - Ai(s) = ^sec^(6')w3(s), 



K2(s) - Aafs) 



^ sec(^) tan(6')a;3(s) H sec^(6')c<;4(s) 

4 16 



K3(s) - A3(s) = - ^ sec^(6')wi(s) - ^ sec(6') tan(6')(x;2(s), 

K4(s) - A4(s) = - -^sec^(6')u;2(s). 
lo 



(22) 
(23) 

(24) 
(25) 
(26) 
(27) 



The symmetry of G [namely G{x, s) = &{s, x), where f denotes the complex conjugate] 
yields a further condition; 

4 4 

^Xj{s)ujj{x) = ^Kj{x)u!j{s). (28) 
i=i i=i 

Because G{x, s) must also be a solution to the adjoint problem G{x, 3)711 = ^2S{x — y) 
(where 1-L\ acts to the left), we see that nj and \j must be linear combinations of the Uj. 
Thus we look for 32 real constants, nl and A^ (where j = 1, 2, 3, 4) which appropriately 
define 

4 

^^^{S)=Y,^^^AS). (29) 
4 

5^Aju;,(.). (30) 



Equations 



Ai(s^ 

then become 
k{ - A{ = (5j3^sec^(6'), 



Kj'-) 



i^i- Xi = -Sj2— sec^{6), 
16 



(31) 
(32) 
(33) 
(34) 
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(a) (b) 




Figure 1. (a) shows the Re [g{x)] and (b) shows the Im [g{x)] as defined by equation 
[35] with A — B — 1. The orange and green fines sfiow the contours of the real and 
imaginary parts of g respectively. 



(where 6jk is the Kronecker deha) while equation [28] becomes 

We can also set Aj = A2 = Ag = A4 = since these diagonal elements only affect the final 
solution for ipi^x) by adding a constant times ojj{x) which is of no physical interest since 
it is just transforming the solution into one of the four previously-mentioned symmetry 
groups. This leaves us with 26 equations for the 32 unknowns, the remaining 6 equations 
are provided by the boundary conditions on ipi{x). 

5. Example problem 

5.1. ID Bose-Einstein condensate in the presence of fluctuations 

Thermal and quantum fluctuations in a Bose-Einstein condensate cause a small- 
but-finite population of non-condensed particles. When a soliton is present in the 
system these non-condensed particles bunch up in the low-density region around the 
soliton [39l dO]. Without paying close attention to the specific details of this non- 
condensed density, we assign g{x) [of equation [7] the following fairly generic form; 

g{x) = cos^(6') [Atanh (xc) sech^ (xc) + iB sech^ (xc)] (36) 

where A and B are real constants [g{x) is shown in Fig. [T]with A = B = 1]. Note that 
we have chosen g{x) to have the same symmetry as ipo (that is the real part is odd, 
while the imaginary part is even) and that g{x) decays at the same rate as 1 — I'l/'oP- As 
boundary conditions on ipi we simply say that ipi{x) — )-constant and Dxipi{x) — )■ 
as X — 7- 00, aswell as basic symmetry arguments; Re['?/'i(x)] = — Re ['i/'i(— a;)] and 
Im ['?/'i(x)] = Im ['?/'i(— x)]. 

Divergences in ■i/'i as x — )■ 00 can be avoided by the conditions; 

A^ = —Ki, A2 = —1^2^ ^3 ~ ~^3' l^"^) 
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and the symmetry is ensured by the conditions; 

= 0, \l = sec\e), Xl = -- sec(^) tan(^). (38) 
4 8 

These six additional conditions give us the Greens matrix, 
- , , sec^(^) / x_ / X sec(6') tan(6') , . , . sec^(6) . . . , 

^11(0; > S) = -^UJi{x)uJ3{s) + '-^ '-^C02ix)uJ3{s) + ' ' UJ2{x)uj4{s) - 

4 o dz 

sec^(6') , ^- r \ sec(^) tan(6') , s- , ^ sec^(^) 

^ W3(x)Wi(s) UJ3{X)UJ2{S) ^;^UJi{x)uJ2{s), 

~ . , sec^(6') , X _ , X sec(^)tan(^) , , . sec^(6) . . , 

Gii(x <S) = UJi{x)UJ3{s) U2{X)UJ3{S) ^^U2{x)UJi{s) + 

sec^(9) / x_ / X sec(6') tan(6') , ^- , ^ sec^(9) 

— - — uj3{x)uji{s) H U3{x)uj2{s) H ^^W4(a;)w2(s), 

~ , , sec^(6') , , , , sec(^) tan(6') , , , , sec^(^) , , , , 

Guix > S) = -^UJi{x)u3{s) + '-^ '-^L02ix)uJ3{s) + —^UJ2{x)u4{s) - 

sec^(6') sec(^^) tan(6') sec'^f^^) 

U3{x)Ui{s) U3{X}U2{S) ^^U4{x)U2{s) , 

~ , . sec'^iO) . . , , sec(^^) tan(^^) , , , , sec^(6') , , , , 

Gi2{x <S) = UJi{x}UJ3{s) UJ2{X)UJ3{S) :^—UJ2{x)UJ4{s) + 

sec^(6') . . , , sec(6') tan(^^) , , , , sec^(6') , , , , 

UJ3{x)UJi{s) H UJ3{X)UJ2{S) H ^^W4(X)W2(S), 

G21 and G22 are easily deduced from the symmetry of G. The expression for ipi then 
follows, 

ipi^x) = ^sech2(xc) 2xc(Acos(2^) + fisin(2^)) + sin(0)(2ficos(^) - 

Asin(^))sinh(2a;c) + ^cos{9)[Asm{9) - 2B cos{9)] (39) 

and ipi is plotted in Fig. |31 One can easily check that equation EHl is indeed a solution 
to equation [7| with g{x) defined by equation [361 

6. Conclusion and discussion 

In this article we have introduced four exact analytic solutions to the NLS equation 
linearised around a dark soliton [equation [T2]. These solutions are given in equations [T^ 
[T71 These four solutions provide a possible means of bypassing the need to solve the 
spatial perturbative correction (denoted ipii^x) in this paper) using the complete set 
of finite E eigenf unctions [given in equations [TOl - [TT] supplemented with generalised 
eigenf unctions for the nullspace of Hx, (a procedure which appears to be common-place 
in the previous literature in-spite of it's apparent difficulty [El [15]). To illustrate this 
point, we constructed a Green's matrix which can be used to find a solution to equation [7| 
once boundary conditions have been defined. We applied the technique to the problem 
of thermal and/or quantum fiuctuations within a Bose-Einstein condensate. 

It is interesting to note that, of the four solutions presented in equations [T4l - [T7[ only 
two of them [oJi{x) and ^2(3^)] remain bounded in the limit as a; — )■ 00. The other two. 
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Figure 2. (a) and (b) show the real and imaginary parts of Gn and (c) and (d) show 
the real and imaginary parts of 6*12 (we have set 6 = 7r/4). The orange and green lines 
show the contours of the real and imaginary parts of the function respectively. 




Figure 3. (a) shows the Re [-01(2;)] and (b) shows the Im [-01(2;)] as defined by equation 
[551 with A = B — 1. The orange and green lines show the contours of the real and 
imaginary parts of respectively. 
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u)3{x) and u)4^{x), are linearly diverging and exponentially diverging respectively. This 
then begs the question as to which set of perturbing functions [g{x) in equation [7] 
are amenable to the use of the Greens matrix defined by equation [191 particulary 
when the boundary conditions require ipi to be bounded. Certainly in the example 
problem of Section 15.11 where the perturbing function itself is strongly localised around 
the soliton, satisfying boundary conditions does not seem to be an issue, since the 
integral in equation [20] is able to contain the divergences associated with and U4^. 
It is also possible to contain divergences by exploiting even or odd symmetries of g{x), 
since and U4 have even and odd symmetries in the real and imaginary parts, the 
integration in equation [20] can once again, avoid undesired divergences. Intuitively one 
might expect (due to the fact that the only interesting parts of equations [T314T7I are in 
the region close to the soliton) that any perturbing function which has a considerable 
nonzero component far away from the soliton would require the use of the radiative 
solutions given in equations [TOl - [TT| and one would follow the procedure of Ref . [I5] . 
However, a general theory on this issue is currently lacking. 
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